function VARbs = doProxySVARbootstrap_single_biascorrected_MSV_CS(VAR,nboot,clevel)
% Constructs bias corrected bootstrapped CIs 

clevel1 = clevel(1);
clevel2 = clevel(2);

jj=1;  %jj=bootstrap draw
res = detrend(VAR.res,'constant');
while jj<nboot+1
        
    
        % Resample full cross-section of residuals within a time period.
        % Jointly draw all residuals in a given year.        
         Tmax = max(VAR.Ti);                                             
         cal  = (min(VAR.Year):1:max(VAR.Year) )' ;                       
        rr    = 1-2*(rand(Tmax,1)>0.5);                                      
        resb = zeros(VAR.T,VAR.Ny);
        for nn =1:VAR.T;
            resb(nn,:)    = rr(find(cal==VAR.Year(nn))) .* res(nn,:) ;  % draw residuals
            VARBS.m(nn,:) = rr(find(cal==VAR.Year(nn))) .* VAR.m(nn,:); % draw proxies
        end
        
                      
           t  = 1;
           Yb = zeros(VAR.T,VAR.n);
           Xb = zeros(VAR.T,size(VAR.X,2)+1); 
           
    for i=1:max(VAR.ii);                                                   % loop over each country        
        Xb(t,:) = [VAR.X(t,:) , 1];                                        % first row for each country set to actual value    
            for j=t:(t+VAR.Ti(i)-1);                                       % loop over time within each country
                 Yb(j,:)   = Xb(j,:)*VAR.bet_bc  + resb(j,:);                 
                 if j<t+VAR.Ti(i)-1;
                 Xb(j+1,:) = [Yb(j,:), Xb(j,1:((VAR.p-1)*VAR.n)) , Xb(j,(VAR.p*VAR.n+1):end)];  % FIXED EFFECTS
                 end
            end 
         t = t + VAR.Ti(i);
    end    

   % run Proxy VAR on bootstrapped data to get Proxy IRF
    % transfer information from VAR
    VARBS.I_m= VAR.I_m; VARBS.k=VAR.k; VARBS.p=VAR.p; VARBS.Ny=VAR.Ny; VARBS.N=VAR.N; VARBS.irhor=VAR.irhor;
    VARBS.normalize=VAR.normalize;
    VARBS.Y  = Yb;
    VARBS.X  = Xb(:,1:end-1);                      % end-1 since X defined without constant
   VARBS.bet = [VARBS.X ones(VAR.T,1)] \ VARBS.Y;
   % bias correct betas
   VARBS.bet_bc = [VARBS.bet(1:(VAR.Ny*VAR.p),:) - VAR.bc_diff ; VARBS.bet((VAR.Ny*VAR.p+1):end,:)];
   % correct alphas
     uhat_bs     = Yb - Xb(:,1:VAR.Ny*VAR.p)*VARBS.bet_bc(1:VAR.Ny*VAR.p,:);
     D1       = [VAR.D(:,1:(VAR.N-1)) ones(VAR.T,1)];
     alpha_bc = (D1'*D1)\(D1'*uhat_bs);
   VARBS.bet_bc(VAR.Ny*VAR.p+1:end,:) = alpha_bc;
   VARBS        = doProxySVAR_biascorrected_single_EV(VARBS);                           
     
   for j=1:VAR.k;
       irs = VARBS.irs(:,:,j);
       VARbs.irs(:,jj,j) = irs(:);
   end
   bsRM(:,jj)        = VARBS.RM ;     
   VARbs.bet(:,:,jj) = VARBS.bet;
   jj=jj+1;   
end  
    
 for j=1:VAR.k  
 VARbs.irsH1(:,:,j)=reshape(quantile(VARbs.irs(:,:,j)',(1-clevel1/100)/2),VAR.irhor,size(VARBS.irs,2));
 VARbs.irsL1(:,:,j)=reshape(quantile(VARbs.irs(:,:,j)',1-(1-clevel1/100)/2),VAR.irhor,size(VARBS.irs,2));
 VARbs.irsH2(:,:,j)=reshape(quantile(VARbs.irs(:,:,j)',(1-clevel2/100)/2),VAR.irhor,size(VARBS.irs,2));
 VARbs.irsL2(:,:,j)=reshape(quantile(VARbs.irs(:,:,j)',1-(1-clevel2/100)/2),VAR.irhor,size(VARBS.irs,2));
 VARbs.irs50(:,:,j)=reshape(quantile(VARbs.irs(:,:,j)',.5),VAR.irhor,size(VARBS.irs,2));
 end
VARbs.RMci=[quantile(bsRM',(1-clevel1/100)/2), quantile(bsRM',1-(1-clevel1/100)/2)];

VARbs.Xb = Xb; %example of last draw
end



